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ABSTRACT 

We develop analytic and numerical models of the properties of super-Eddington stellar winds, 
motivated by phases in stellar evolution when super-Eddington energy deposition (via, e.g., 
unstable fusion, wave heating, or a binary companion) heats a region near the stellar surface. 
This appears to occur in luminous blue variables (LBVs), Type Iln supernovae progenitors, 
classical novae, and X-ray bursts. We show that when the wind kinetic power exceeds Edding¬ 
ton, the photons are trapped and behave like a fluid. Convection does not play a signihcant role 
in the wind energy transport. The wind properties depend on the ratio of a characteristic speed 
in the problem VcUt ~ (£G)^^^ (where E is the heating rate) to the stellar escape speed near 
the heating region Vsscirh)- For VcUt ^ i'esc(?‘/i) the wind kinetic power at large radii E„ ~ E. 
For Vcrit ^ Vesc(i‘/i), most of the energy is used to unbind the wind material and thus E„ < E. 
Multidimensional hydrodynamic simulations without radiation diffusion using ELASH and 
one-dimensional hydrodynamic simulations with radiation diffusion using MESA are in good 
agreement with the analytic predictions. The photon luminosity from the wind is itself super- 
Eddington but in many cases the photon luminosity is likely dominated by ‘internal shocks’ 
in the wind. We discuss the application of our models to eruptive mass loss from massive stars 
and argue that the wind models described here can account for the broad properties of LBV 
outflows and the enhanced mass loss in the years prior to Type Iln core-collapse supemovae. 
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1 INTRODUCTION 


During most phases of stellar evolution, the structure of stars can 
stably adjust so that energy generation in the stellar interior is less 
than or of order the Eddington luminosity, ensuring that the star is 
on average in thermal and hydrostatic equilibrium. However, this 
balance can be upset by instabilities or external perturbations (e.g., 
a binary) leading to epochs of super-Eddington energy generation. 
The canonical example of this phenomena is runaway thermonu¬ 
clear fusion due to either fusion under degenerate conditions (e.g 
Mestel lj52|l or gas pressure dominated thin shell fusion (e.g 


Schwarzschild & Harm|1963)). This occurs e.g., during the He shell 


flash in low mass stars. He shell fusion on the Asymptotic Giant 
Branch, classical novae, radius expansion X-ray bursts on accret¬ 
ing neutron stars, and Type la supemovae. 


In contrast to the case of mnaway fusion, however, the prop¬ 
erties of super-Eddington heating at a given location in a star need 
not be solely determined by the stellar conditions at that location. 
Instead, energy can effectively be deposited at a given radius by an 
external source (e.g., via tidal heating or via a companion during 
common envelope evolution; e.g., |Paczynski|197^ or by non-local 
redistribution of energy (e.g., via wave transport of energy in stellar 


interiors; e.g., |Piro|201 l[|Quataert & Shiode|201^ . Eorthis reason, 
we shall largely use the terms ‘energy generation’ and ‘energy de¬ 
position’ synonymously in this paper. 


There is strong evidence that massive stars undergo periods of 
super-Eddington energy generation/deposition, although the phys¬ 
ical causes are much less well understood. Luminous blue vari¬ 
ables (LB Vs) such as Eta Carinae radiate a photon luminosity sig¬ 
nificantly exceeding the Eddington luminosity for months-decades 
(many dynamical times) and drive a wind whose time-averaged ki¬ 
netic power exceeds both the Eddington luminosity and probably 
the photon luminosity jSmith et al.||2003l see, e.g., [Davidson &| 
|Hum phre ys|2012| for a review of LBVs). Such outbursts may dom¬ 
inate the total mass loss from massive stars (e.g., [Smith & Owockij 
[2006[|Kochanek[201 1[ ). Moreover, ~ 10% of supernova (SN) pro¬ 
genitors experience enhanced mass loss in the decades to weeks 
prior to core collapse (much larger than can be explained by line 
driven winds). Evidence for this powerful mass loss includes obser- 


vations of luminous outbursts that precede supernovae dPastorello 

et al. 

|2007[ [Foley et al.[|2007[ [Mauerhan et al.||2013[ [Ofek et al. 

2013 

and mass-loss rates ~ lO^^* - 1 Mq yr“‘ inferred from obser- 
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vations of circumstellar interaction in Type Iln SNe (e.g., |Kiewe| 
|etal.|2012|[Smith|2014| ). 

What is the response of a star to super-Eddington energy gen¬ 
eration/deposition? There are essentially three regimes. If sustained 
energy deposition leads to heating on a timescale short compared 
to the local dynamical time, a propagating shock will form jDessart| 
|et al.||2010| ). By contrast, if the energy deposition occurs on a 
timescale long compared to the dynamical time, the response is 
at least initially largely hydrostatic: the increase in thermal energy 
slowly lifts matter out to larger radii and convection sets in in re¬ 
sponse to photons being unable to carry the energy (e.g., [Joss et al.| 
|1973| >. Finally, if the super-Eddington energy generation occurs for 
sufficiently long and/or sufficiently close to the stellar photosphere, 
it can drive a powerful wind (e.g., |Kato & Hachisu|1994[|Owocki| 
|et al.|200^ . In this paper, we are interested in the latter regime. This 
requires that the energy deposition is on for several thermal times 
so that the entire star exterior to the heating region inevitably ad¬ 
justs its structure significantly in response to the energy deposition. 
It is important to stress that this is not always the case. For example, 
during the Helium flash in low mass stars, runaway fusion produces 
a highly super-Eddington energy generation rate that locally drives 
vigorous convection. However, the expansion of the exterior of the 
star in response to the energy input from fusion lifts the degeneracy, 
quenching the fusion. There is no noticeable increase in the surface 
luminosity (it in fact decreases slightly; e.g., |Bildsten et al.|2012| ), 
let alone a powerful wind. 

Previous work on super-Eddington stellar winds has in many 
cases approached the problem as a radiative transfer problem. For 
example, increases in opacity at specific temperatures (e.g., the iron 
opacity bump) can produce a locally super-Eddington flux and po¬ 
tentially drive a wind (|Kato & Hachisu|1994 JEichlCTe^^alJl^^ 


|Heger & Langer 1996 1 . Alternatively, |Shaviv| i 2001^ ; |Owocki et ak] 




van Mar 


et al.|p008ll argued that super-Eddington winds 


are regulated by how radiation dilfuses through a highly porous 
stellar atmosphere that is envisioned to be the outcome of some 
combination of convection, (magneto)-hydrodynamic instabilities, 
and/or the wind launching process itself. In this paper, we show that 
when the wind kinetic power exceeds the Eddington luminosity, the 
photon diffusion time is long compared to the advection time in the 
wind (^3^: photons are trapped and behave like a fluid. As a result, 
we argue that how the radiation diffuses through the wind/stellar at¬ 
mosphere is less dynamically important than suggested in previous 
calculations. Instead, super-Eddington winds are essentially hydro- 
dynamic, more analogous to |Parker| ( |1958^ ’s thermal solar wind 
model (generalized to a radiation pressure dominated fluid) than 
standard radiation-driven winds from massive stars. 


1.1 Outline of This Paper 

We study the response of a star to continuous energy deposition 
near the surface at a rate larger than the typical Eddington lumi¬ 
nosity. Throughout much of the paper, we are relatively agnostic as 
to the physical origin of this heating. It could represent, e.g., the 
outcome of unstable nuclear fusion or heating via waves generated 
in the core that carry a large energy flux to the surface where the 
waves dissipate or heating triggered by a binary companion (e.g., 
tides, common envelope). 

As noted above, the response of a star to super-Eddington en¬ 
ergy deposition has been studied extensively in the context of run¬ 
away fusion in stars. This work has tended to focus on the limit in 
which sustained heating for of order a thermal time or more leads 
to the formation of an extended convection zone that carries the ex¬ 
cess energy to larger radii (e.g.,|Woosley et al.|2004[|Weinberg et al.| 


|2006[[Prro & Chang|2008) . We briefly review this quasi-hydrostatic 
response in ^ We then discuss the role of convective energy trans¬ 
port in wind solutions, rather than hydrostatic envelopes ( §3.1^ . We 
show that energy transport by convection is subdominant in steady 
state winds. Thus although convection can play an important role 
inflating a stellar envelope in response to energy deposition, it be¬ 
comes unimportant once a roughly steady wind develops. This typ¬ 
ically occurs on a few thermal times. 

In ^3.2| we introduce and solve a simple model problem for 
super-Eddington winds neglecting energy transport by convection 
and radiation. These are radiation driven winds, in a regime where 
the photons are trapped and behave like a fluid, ^presents numer¬ 
ical examples of wind generation using both FLASH and MESA 
that are in good agreement with the analytic models, ^presents 
some of the observational properties of super-Eddington stellar 
winds. Finally, discusses several applications of our work to 
massive stars and 0 summarizes our main conclusions and high¬ 
lights some key questions not addressed by our models. 

A central feature of all of the models in this paper is that we 
assume that an unspecified process (e.g., fusion, wave dissipation) 
leads to heating at a rate E near a heating radius radius rn. This 
heating is put in ‘by hand’ in our calculations. The advantage of 
this treatment is that it means that our model is potentially appli¬ 
cable to a range of physical processes and stellar contexts. How¬ 
ever, when scaling our analytic and numerical results, we focus on 
the application to super-Eddington winds from massive stars. It is 
also important to stress that the photon luminosity is not an input 
quantity in our model, as is often the case when super-Eddington 
wind models are formulated in terms of what wind properties re¬ 
sult from a given super-Eddington luminosity (or, equivalently, a 
given Eddington ratio T). Instead, the photon luminosity is a de¬ 
rived quantity in our model, which depends on the heating rate E 
and the properties of the resulting wind (^. We believe that this is 
a more appropriate formulation of the problem of super-Eddington 
stellar winds for the systems of interest in this paper. 

2 HYDROSTATIC ADJUSTMENT IN RESPONSE TO 
ENERGY DEPOSITION 

Assume that there is some source of heating at a rate E near a heat¬ 
ing radius rh.The local thermal energy increases in response to the 
heating on the timescale ~ Anr^uHlE where u and H are the ther¬ 
mal energy per unit volume and scale-height at rj,, respectively. We 
assume throughout that the source of heating is on for multiple ther¬ 
mal times so that the outer envelope of the star inevitably adjusts 
its structure significantly in response to the energy deposition. 

After of order one local thermal time, the stellar envelope be¬ 
gins to expand outwards due to the increased thermal pressure. 
Neglecting for now the possibility of a wind, super-Eddington en¬ 
ergy deposition also inevitably drives convection, which attempts 
to carry the energy that photons cannot (e.g., |Joss et al.|1973| >. In 
hydrostatic models (i.e., absent winds), this leads to the creation of 
a large convective envelope even if the star was initially compact. 
The total timescale to rearrange the structure of the stellar envelope 
is of order the thermal time of all of the mass exterior to ~ rj,: 

GM{<r^)M{>n) 

Hhermal “ 

2r^E 

where M{< ri,) is the stellar mass within rj, and M{> r^) is the stel¬ 
lar mass exterior to rn. The thermal time as defined in equation^ 
is often significantly longer than the local thermal time at ri, be¬ 
cause all of the stellar envelope participates in the convective and 
hydrostatic readjustment. 
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Figure 1. Density profile as a function of radius at different times after the 
onset of energy deposition with £ = 3x 10® Lq at r;, ^ 2Rq (in an 11.1 Mq, 
Z = Zq model at core He exhaustion). These are hydrostatic MESA models 
that cannot develop a wind. The initial thermal time (eq. ^ at r;, is ^ 0.2 
yr, which sets the initial expansion time of the envelope. At late times the 
heating has generated a spatially extended convective envelope with a den¬ 
sity profile similar to the p oc profile expected for radiation dominated 
convection. Compare with Figure|^which shows analogous density profiles 
in hydrodynamic MESA models that do develop a wind. 

To quantitatively illustrate this process, we have carried out 
a number of calculations using the MESA stellar evolution code 
( [Paxton et al.|2011||2013|[20T5l > in which we inject energy into the 
star at a specified rate. The inlists for the models in this section 
are given in Appendix jAlj The two key properties of these models 
are, first, that they are essentially hydrostatic, so that there is no 
option for a wind to develop. Secondly, for simplicity we utilize 
the MLT-l-l- option in MESA which forces the convection to he 
efficient, i.e., to have roughly constant entropy, in radiation pressure 
dominated regions ( [Paxton et al.|2013] l. This can formally become a 
poor assumption at large radii, but we show in the next section that 
once a wind develops convection in fact plays little role in setting 
the wind properties. 

We focus here on an initially compact blue supergiant (BSG) 
because this highlights most dramatically how the stellar structure 
adjusts in response to energy deposition in the envelope. The stellar 
model is derived from a 30 Mq, Z = Zq star evolved to He core 
exhaustion without energy deposition. At He core exhaustion, the 
stellar mass and radius are 11. IMq and 2.1 Rq, respectively. We then 
deposit energy at a rate E = 3 x 10®Lq in a region centered at 
r\i = 2i?o- This choice of E is approximately 10 times the electron¬ 
scattering Eddington luminosity. The stellar mass exterior to th is 
initially M(> r*) ~ 3 x IO^^Mq although this increases to M{> 
r;,) ~ 0.03Mo during the hydrostatic readjustment. 

Figure^ shows the density profile of the star at several differ¬ 
ent times as it expands from the initially compact BSG to become 
a red supergiant (RSG). The expansion occurs on the thermal time 
in equation^ this is 0.2 yr for the initial BSG model but the en¬ 
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Figure 2. Convective (Lconv) and total (Ltot, sum of radiative and convec¬ 
tive) luminosities as a function of radius at dilferent times after the onset 
of energy deposition with E = 3 x lO^L© at r/, ^ 2Rq (in an 11.1 Mq, 
Z = Zq model at core He exhaustion). These are hydrostatic MESA models 
that cannot develop a wind. The initial thermal time (eq. ^ at r/, is ^ 0.2 
yr, which sets the initial expansion time of the envelope. At late times most 
of the energy is carried to lai'ge radii by convection. L„jo.v = 4;rr^pc‘'] is the 
maximal convective power for subsonic convection (not applicable in these 
models because of the use MLT++ in MESA; see ©■ Compare with Eig- 
ure [10[ which shows analogous luminosity profiles in hydrodynamic MESA 
models that do develop a wind. 

tire process of inflating the stellar envelope takes about a factor of 
~ 10 longer as more mass expands out and is incorporated into the 
convective envelope. The density profile at large radii in the RSG 
configuration is well described by a power-law. This is very close to 
a y = 4/3 polytropic atmosphere expected for efficient convection 
in a radiation pressure dominated fluid (as we derive analytically in 
the next section). 

Figure]^ shows the total (radiative and convective) Ltot (black) 
and convective Lconv (red) luminosities as a function of radius for 
the same models shown in Figure[T] For comparison, we also show 
the energy deposition rate of £ = 3 x 1O®L0. Recall that the en¬ 
ergy deposition occurs at rh - 2 Rq. This is why L,ot and Lco„v rise 
abruptly at that radius. Figure|^shows that at early times (1 = 0.2 
yr), very little of the energy deposited into the star has gone into 
photon or convective power (since Ltot.Tconv « E). Instead, most 
of the energy has gone into heat and P-dV work to expand the en¬ 
velope. At later times, however, the star has reached a new equilib¬ 
rium in which much of the energy supplied at small radii is carried 
to large radii via convection and, to a lesser extent, photon diffu¬ 
sion. 

3 ANALYTIC MODELS OF SUPER-EDDINGTON WINDS 

In this section we derive the properties of steady state spherically 
symmetric radiation-pressure driven winds in response to super- 
Eddington energy deposition. We first show that convection is not 
an important energy transport mechanism in such winds, in spite 
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of its role in initially expanding the stellar envelope in response 
to energy deposition (as shown in We solve a model steady 
state spherically symmetric wind problem in §3.2| and describe 
how the wind properties depend on the energy deposition rate E 
and the properties of the stellar model. Throughout we neglect dif¬ 
fusive transport of energy and assume that the optical depths are 
sufficiently high that the photons are trapped and so behave like 
ay = 4/3 fluid. Physically, this assumption is motivated by the 
fact that photon diffusion can only transporf energy at the Edding¬ 
ton luminosity and yet we are interested in problems for which the 
energy deposition rate is super-Eddington. We check the validity 
of this assumption in ^3.3| and show that it requires that the wind 
kinetic energy flux exceed the Eddington luminosity. In ^3.4| we 
briefly contrast our models with super-Eddington wind models for¬ 
mulated in terms of the Eddington ratio E. 


3.1 The Role of Convection in Wind Solutions 


The blue line in Eigure|^shows the maximum power that subsonic 
convection can carry for the hydrostatic MESA model at t = 3.4 
yr. The convective energy transport should physically be limited by 
the fact that the velocities remain subsonic, so that 


L 


max 


~ And pel. 


( 2 ) 


Eigure shows that the maximum convective power decreases 
rapidly at large radii and by r ~ 50Rq convection would have diffi¬ 
culty carrying the required energy flux (this is not true in the numer¬ 
ical models because of MLT++ in MESA). Assuming the convec¬ 
tive power Lconv ~ E is super-Eddington, photons also cannot carry 
the energy outwards. An a priori plausible hypothesis is that be¬ 
cause of the failure of either convection or photons to carry the en¬ 
ergy E outwards, the pressure will build up at roughly the location 
where convection ceases to efficiently carry the energy outwards. 
The radiation pressure gradient would then accelerate the matter 
outwards in a wind. We now show, however, that this hypothesized 
transition cannot occur in a steady state wind and instead convec¬ 
tive energy transport is sub-dominant throughout the wind. 

Consider a hypothetical wind in which the flow is heated near 
a radius rj, at a rate per unit volume q (with f dr And q - E) and 
then driven outwards by radiation pressure at radii much greater 
than the heating radius rn because of a failure of convection to 
transport the energy outwards. Neglecting radiation transport of en¬ 
ergy, the energy equation describing such a flow in steady state and 
spherical symmetry would be 


■ V7 f ■ ^ ^Eonv -T., 

pvT— = q-V ■ = q -—. (3) 

dr Anr-^ dr 


response to energy deposition (see convection is unimportant 
once a steady wind is established. For these reasons, we do not con¬ 
sider convective energy transport in our wind solutions that follow. 
The numerical models in ^confirm the sub-dominance of convec¬ 
tion in the wind dynamics. 


3.2 Spherically Symmetric Super-Eddington Winds 


We consider a simplified model problem to describe the physics of 
stellar winds with super-Eddington energy deposition. The struc¬ 
ture of the stellar envelope can be modeled analytically if we ne¬ 
glect its self-gravity. We take P = . where is an ef¬ 

fective entropy of the envelope and 7 c„y quantifies its stratification. 
Solving hydrostatic equilibrium, dPfdr = —GpMjr^, where M is 
the total stellar mass interior to the stellar envelope, yields the den¬ 
sity and sound speed profiles: 


Penv(r) — 


1 \l/(7env A) 

yenv 1 \ / 

yenv^env/ ' ^ 


~R~ ) 


(4) 


and 



yenv - 1 /GM _ GM\ 
Tenv V r R I 


(5) 


The stellar radius R in equations|^&|^is defined fo be where p = 0. 
Note also that in equationj^the sound speed c, is defined without a 
factor of the adiabatic index, so it is more analogous to the isother¬ 
mal sound speed. We use this definition throughout the paper. 

The analytic density profile in equation|^for a radiation domi¬ 
nated convective envelope with yenv = 4/3 becomes p oc r^^, which 
is in good agreement with the numerical results in Figure ^ for 
1 = 0.9 and 3.4 yr, i.e., once R » r\,. The small difference in 
power law slope is because in the numerical model gas pressure 
contributes about 10% of the total pressure, which modifies the adi¬ 
abatic index from the pure radiation pressure value. 

Given equation|^as the background density profile of the stel¬ 
lar envelope, we solve the steady state spherically symmetric wind 
equations for a wind subject to a total heating rate E at radius r*. In 
the analytic calculation it is convenient to assume that the width of 
the heating region is « p, so that the heating is relatively spatially 
localized. Mass, momentum, and energy conservation for a wind in 
thermal equilibrium (equal radiation and gas temperatures) can be 
written as 


M = 47rr“pv = constant 


( 6 ) 


dv _ \dP GM 
dr p dr r- 


(7) 


At large radii » q,, ^ = 0 by assumption. If the failure of 
convection to transport energy outwards is to drive a wind then 
dLcomIdr < 0 as energy is transferred from convection to thermal 
energy and then to a wind. This implies, however, that the right- 
hand side of equationj^is > 0 and thus that dsfdr > 0, i.e., the flow 
is convecfively stable. This demonstrates that there is a not a consis¬ 
tent steady state solution in which the flow transitions with radius 
from a large hydrostatic convective envelope to a wind. Moreover, 
once a steady wind is established, convection is not an important 
source of energy transport. To see the latter, again consider equation 
[^but now with Fconv = 0- In the heating region ^ > 0, ds/dr > 0 
and so the flow is convecfively sfable. Oufside fhe healing region, 
q = 0 and the flow is buoyantly neutral. 

The above considerations demonstrate that although convec¬ 
tion does play a role in initially expanding the stellar envelope in 




a , gmw 

(f-rad + M Bej - — 

End + M 

-v^ + h - 

(2 r j 


where is the total power carried by photons. Be in equation]^ 
is the Bernoulli parameter, the conserved energy per unit mass flux 
for a steady wind absent radiation, and h is the enthalpy of the fluid. 

To solve equations |6|8| analytically, we make several simplify¬ 
ing approximations. First, we neglect because we are interested 
in systems with super-Eddington energy deposition. We assess the 
validity of this approximation in §3.3| We further assume that out¬ 
side the heating radius, the fluid is radiation dominated with an 
adiabatic index y = 4/3. This is appropriate even if the stellar en¬ 
velope prior to heating has a stratification ye,„ t 4/3. The reason 
is that if the net heating in the heating region is super-Eddington 
it will increase the entropy to the point where radiation pressure 
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dominates. For a fluid with adiabatic index y = 4/3, the enthalpy 
in equation]^ is given hy h = 4Cj, where = Pjp and P is the 
radiation pressure. 

Equations 1^ and are equivalent formulations of the steady 
state, spherically symmetric energy equation, but the latter is more 
convenient for the wind problem. The reason is that equationj^can 
be integrated to yield global conservation of energy in the wind 

M [Be{r^) - Be(ro)] = M [Be{r oo) - Be(ro)] = E (9) 


where Be{rk) is assumed to be evaluated just outside the heating 
region and is a radius in the stellar envelope just interior to 
where the heating occurs. The first two expressions in equation 
are equivalent because q is assumed to be zero outside ~ r* and 
hence energy is conserved between r;, and large radii. Beir^) in 
equation]^ is effectively the binding energy per unit mass of the 
stellar envelope prior to heating. We define 


Be(ro) = ( 10 ) 

Physically, equation[TO|corresponds to the assumption that a typical 
speed of order the escape speed at the heating radius is required to 
unbind matter from the stellar envelope. The dimensionless param¬ 
eter / quantifies the binding energy of matter in the stellar envelope 
and depends on the exact structure of the envelope. For example, for 
a polytropic atmosphere with negligible mass (so that self-gravity 
can be neglected) described by equations |^ & |^ it is straightfor¬ 
ward to show that Be(r„) = -GMIR for yem = y (independent of 
both radius Vo and yenv)- Thus for stellar envelopes with negligible 
mass it is possible to have / <K 1 if » r;,. For more realistic 
massive stellar models, we find that / ~ 0.1 - 1 in the outer stel¬ 
lar envelope, with smaller values for more extended envelopes (see 
Fig .|1 l|discussed in 

Equations |6|8| can be combined to yield the wind equation that 
describes the acceleration of the wind (again neglecting Lrad): 


I dv / 2 4 2 \ 8 q GM 
vdr\ 3 '7 3 r 3vp P 


( 11 ) 


Equation [m is of the usual form where the requirement that the 
flow smoothly evolve from subsonic to supersonic at a sonic point 
implies that the left and right hand sides of equation [TT] vanish so 
thaQ 


Sonic Point —> 


7 

V“ 


and 


8c? 

3 c. 


g 

3vp 



( 12 ) 


Equations [T^ are additional boundary conditions that specify the 
mass outflow rate M and the location of the sonic point r,. The solu¬ 
tion of equation[T^depends on the relative value of two timescales 
at the sonic point, the heating timescale 4eat ~ pdjq relative to 
the dynamical timescale fdyn ~ rjcs- We are working in the limit 
fheat <: fdyn> Of else the heating would be dynamical and drive 
shocks. In this limit the solution of equation[^for the sonic point 
is the usual sonic point condition absent heating: 




3 GM 
8 r. 


(13) 


To solve for the wind mass loss rate M and sonic point radius 
Ts, we begin by noting that because the sonic point is located at 


* There are also breeze solutions that never go super-sonic. These may be 
relevant to the early time dynamics in the numerical models presented in 
^ when the envelope mass just exterior to the heating region is so large 
that it stifles the wind. 



Figure 3. Analytic estimates of the properties of steady state spherically 
symmetric super-E ddin gton winds (see asymptotic wind speedva, (in 
units of Vcrit; see eq. |22^ , mass loss rate M (in units of Merit: see eg. [23} , and 
asymptotic wind energy flux (in units of the energy input rate E). 


rs> Vf,, energy conservation implies fie(rj) Be{r —> co) v?„/2, 
where v„ is the asymptotic velocity of the flow at large radii. Using 
equation[T^ Beir,) can be shown to be 3GMI4rs, which implies 


3 GM 


(14) 


An expression for the asymptotic velocity can be obtained by 
combining equations [^&|10| 

vi = ^ - fvi,(rp. (15) 

M 

The mass loss rate can be evaluated at the sonic point as 


M = 4;rr?p(f'j)v(rj) = dddn 


{GMfpir,) 


(16) 


where we have used equations|12|and[T4| To simplify equation [T^ 
further, we require the density at the sonic point pfr,). Since the 
flow is subsonic interior to the sonic point, this can be estimated to 
reasonable accuracy by solving hydrostatic equilibrium interior to 
the sonic point with a boundary condition at small radii that deter¬ 
mines the normalization of the density. To do so, we proceed as fol¬ 
lows. Eor radii satisfying rj, < r < there is no heating and so the 
outflow is adiabatic with y = 4/3 by assumption. Moreover the flow 
is subsonic and so roughly in hydrostatic equilibrium. As a result, 
the density profile between the heating region and the sonic point 
is well described by equation|^with yem 4/3, and 

» oo (the latter because the wind is unbound). This implies 

1 (GMP 

(17) 

'"wind 


The wind entropy /twind can be calculated by using pressure bal¬ 
ance across the heating radius. We define and to be radii just 
inside and outside the heating radius, respectively. Pressure bal¬ 
ance implies P{r,) r. P(r+) = ffwmdP(r+)"/' = (GM^ 
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Table 1. List of wind models evolved with FLASH, and summary of results. Columns (left to right) show model name, heating rate, initial envelope radius, 
radial width of the heating region, dimensionality, ratio Vesc(^)/vcrit (eq. |22^ , total simulated time, time to achieve steady-state in the sonic point evolution 
(dlnr,./dlnt = 10“^), steady-state mass loss rate in absolute units and normalized to the critical value (eq. |23t , ratio of steady-state outgoing wind energy 
loss rate to input power, and ratio of the asymptotic sonic point radius to the initial envelope radius. For uniformity, the mass and energy loss in 2D and 3D 
have been multiplied by 4;r/An, where AH is the solid angle subtended by the computational domain. Results have been rounded to two significant digits. 


Model 

E 

(n^PhVtsciri,?) 

R 

(rh) 

An, 

(rh) 

Dim. 

l'esc(^)/^'crit 

^h\m ^^steady 

irh/Ve&ciEh)) 

M 

(rh^PhVescirh)) 

M/Mcnt 

£w/£ 

rJR 

L005R2.5-ld 

5.2E-4 

2.5 

0.075 

ID 

5.4 

2 . 8 E +4 

1.4E+4 

2.5E-3 

3.3E-2 

1.3E-2 

50 

L015R2.5-ld 

1.5E-3 




4.3 

8 . 8 E +3 

5.2E-h3 

7.4E-3 

5.1E-2 

2.9E-2 

25 

L046R2.5-ld 

4.6E-3 




3.5 

5.3E+3 

1.6E-h3 

2.2E-2 

7.9E-2 

6.1E-2 

12 

L150R2.5-ld 

1.5E-2 




2.7 


4.1E-h2 

6.8E-2 

1.2E-1 

1.3E-1 

5.0 

L150R5.0-ld 


5 



2.2 


6.1E-t2 

1.2E-1 

1.7E-1 

1.8E-1 

3.3 

L150R10-ld 


10 



1.6 


5.5E-t2 

2.0E-1 

2.4E-1 

3.1E-1 

1.5 

L150R25-ld 


25 



1.1 


4.1E-h2 

2.9E-1 

3.0E-1 

5.2E-1 

5.4E-1 

L150R50-ld 


50 



8 . 1 E -1 


2.3E-h2 

3.3E-1 

3.2E-1 

6.4E-1 

2.5E-1 

L150R100-ld 


100 



6 .OE -1 


1.8E-H1 

3.6E-1 

3.3E-1 

7.2E-1 

1.2E-1 

L005R2.5-2d 

5.2E-4 

2.5 

0.075 

2D 

5.4 

2.8E+4 

1.4E+4 

2.6E-3 

3.4E-2 

1.4E-2 

51 

L015R2.5-2d 

1.5E-3 




4.3 

8.8E+3 

3.8E-h3 

7.5E-3 

5.2E-2 

3.0E-2 

24 

L046R2.5-2d 

4.6E-3 




3.5 

5.3E+3 

1.6E-H3 

2.2E-2 

7.7E-2 

6.2E-2 

12 

L150R2.5-2d 

1.5E-2 




2.7 


4.0E-h2 

6.7E-2 

1.2E-1 

1.3E-1 

5.0 

L150R5.0-2d 


5 



2.2 


6.1E+2 

1.2E-1 

1.7E-1 

1.9E-1 

3.2 

L150R10-2d 


10 



1.6 


5.5E-t2 

2.0E-1 

2.4E-1 

3.1E-1 

1.5 

L150R25-2d 


25 



1.1 


4.0E-H2 

2.9E-1 

3.0E-1 

5.3E-1 

5.3E-1 

L150R50-2d 


50 



8 . 1 E -1 


2.3E-h2 

3.4E-1 

3.3E-1 

6.4E-1 

2.5E-1 

L150R100-2d 


100 



6 .OE -1 


1.4E-t2 

3.6E-1 

3.3E-1 

7.2E-1 

1.2E-1 

L046R2.5-3d 

4.6E-3 

2.5 

0.075 

3D 

3.5 

2.8E+3 

1.8E+3 

2.2E-2 

7.8E-2 

6.1E-2 

12 


Hydrostatic equilibrium in the stellar envelope requires P(rJ^) - 
pcm{rh)GMHIrl where H is the local pressure scale height in the 
envelope near r^. We define the hydrostatic envelope mass as 

Me„v = WxAn d H (18) 


The factor of 10 in equation [T^ is arbitrary but is motivated by 
the fact that for a yenv = 4/3 envelope H ~ rjA, p^nv and 

thus each decade in radius contributes comparably to the total mass. 
For a spatially extended envelope the total mass is thus a multiple 
~ 10 of the local mass at r;,. Despite this particular motivation, 
equation [T^ can simply be viewed as a definition of used to 
set the density scale for the envelope mass and thus the outflow. 
Note also that if the size of the heating region is comparable to or 
larger than the local scale-height H then P{r^) - P(d ) is not a good 
approximation. Since Ffr/) oc this can be roughly captured 
in what follows by a reduction in the envelope mass M^m- 

With equation^] pressure balance across the heating radius 
implies 


V'l ^ 


Wn (GMf 
43 M 


(19) 


Combining equations^]^]and^]we find the wind density at the 
sonic point of p(r^) == Menv/lOtrrJ and thus arrive at our final 
expression for the mass outflow rate 


M ^ 


4 Menv vj 
15V3 M G 


4 4Te„v 
15 V3 GM 




( 20 ) 


We reiterate that Vesc(r;,) here is the escape velocity of the stellar en¬ 
velope just interior to the heating radius. This is important because 
it sets the Bernoulli parameter of the stellar envelope and thus the 
energy that must be supplied to unbind mass from the star. 

The solution of equation |2Q| depends on the relative magni¬ 
tude of the escape velocity Vcxio,) to a characteristic velocity of 


the problem Vcrit, where 

Vescin) - 620 km s^' 
M = 30 M 30 M0, r* = 30 r;,3o Rq, and 


Vent = 


M 

44env 


1/5 

GE\ == 190km s^ 


M 


103 Me, 


1/5 


71/5 


( 21 ) 


( 22 ) 


where £7 = £'/( 10 ^Lo); uote that 10 ^Lq corresponds to about 10 
times the electron-scattering Eddington luminosity for our fidu¬ 
cial 30 Mq star. In equation]^ and what follows we take a typical 
Menv ~ As discussed in ®(Fig.[T 3 , this is required in or¬ 

der to explain the timescale of observed super-Eddington mass loss 
from massive stars. 

We also define a characteristic mass loss rate using 


E = 


2 4ifcritVe,i, 


Merit — 


2vLt^. 

GM 


(23) 


in which case equation|^can be written in dimensionless form as 


(Merit/ 675 \ /Werit ] 


(24) 


Eigure shows the numerical solution of equations and 
l^for the wind asymptotic velocity Voo in units of Verit, the mass 
outflow rate in units of Man and the asymptotic wind kinetic energy 
flux E„ = 0.5Mv^ in units of the energy input rate E. These results 
show that there are two regimes with dilferent wind physics: 


Regime 1 : Ved, > fF3vesc(rh) 


This corresponds to relatively large stellar progenitors and/or low 
mass, weakly bound envelopes. In this case EjM » f'^lscGh) and 
and|20 is M 0.36 Merit, ke., 


the solution to equations 15 

, 2/5 


M 


0J2(^) £:3/5^iM0yr-‘ 


M 


103 Me, 


- 2/5 


7,3/5 


(25) 
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and 

v..L7«.3<K,kms-'(^)‘'’f;» (26) 

In this regime the sonic point is located at 

/ x-2/5 

r, 0.52 - &0R^ M 30 j ^ 7 "'' 

(27) 

and the asymptotic wind kinetic power is given by 

E. = - E. (28) 

Equation|28|corresponds to nearly all of the energy deposited in the 
stellar envelope going into driving a wind. 

Note that equation|^for the sonic point location can also be 
rewritten as 


r, - 


1 

4 



(29) 


Physically, solutions must have r, > r* since the sonic point can¬ 
not lie interior to the heating radius. This implies that we require 
i'csc(^a) <: 2 vcrit and thus / < 1 for physical wind solutions to both 
be in the regime Vcnt <: E^^Vesdri,) and have > r*. Note that 
Vescdh) <: 2V(,rit and equationtogether imply that the maximal 
value of the asymptotic wind speed is ~ Vesc(r*). 

What happens if Vcru ^ O.SVescfr*) (and hence eg. |29| implies 
Ti ^ ^h)^ This is certainly physically realizable. A simple estimate 
shows that ?ti^ermal(^/7)/^dyn(Oj) (t'esc(^/j)/t'crit) SO that Vcrit ^ ^escdh) 
corresponds to heating on less than or of order a dynamical time, 
which will lead to a strong shock rather than a steady wind. 


Regime 2 : f‘^^Vesc(rh) > Vcnt 


This corresponds to relatively compact stellar progenitors and/or 
higher mass, tightly bound envelopes. In this case the solution to 
equation 20 satisfies EjM » and 


E - {mvL(0 .) ^ ^ - 0.3 Mo yr-' £, M/J r„,3o T' 


and 


^ 1.9 


Erh 


1/3 


200 km s 


-1 / ETri,^3o\ ^ / 

' /Af3o ] 


M 


(30) 


1/3 


In this regime the sonic point is located at 
0.14r;,/^/^(vesc(f^,)/vcrit)‘“^^ i.e.. 


r, ^ 200 Ro£: 


-2/3„5/3 -2/3^2/3 
'■^30 'li30J 


M 




-2/3 


and the asymptotic wind energy flux is given by 


fv^Ari,y 




5/3 


M 


I 103 Me, 


2/3 


(31) 


(32) 


(33) 


In this regime the asymptotic energy flux of the wind is small com¬ 
pared to the energy supplied to the stellar envelope. Most of the 
input energy is used to unbind the gas from the potential of the star, 
as implied by £ \Be(r„)\M O.SfMvl^d^i,) in equation ; 


30 


3.3 Validity of the Adiabatic Approximation 

In deriving the wind properties in this section we have neglected 
photon transport of energy and assumed that the outflow is adia¬ 
batic. This requires that t^,fs ^ dKpjc > fexp - r/v at least out to 


the sonic point {k here is the opacity). Using equations [T^ & 

|16| it is straightforward to show that the condition for adiabaticity at 
the sonic point takes the intuitive form £„, > (this also follows 
directly from equationj^by noting that MBe = £„,). Physically, this 
condition states that if the asymptotic kinetic power of the wind is 
greater than the Eddington luminosity, the photons are trapped in 
the wind at the radii where the wind properties are set (between 
the heating region and the sonic point). The wind properties can 
thus be calculated neglecting radiation diffusion and assuming a 
7 = 4/3 radiation dominated fluid. Since E^jE declines with in¬ 
creasing Vesc(r;,)/Vcrit (Eig.[^& eq.|33[l, in practice the current solu¬ 
tions formally apply only if VKc(r*)/Vcri, < l.lf-^^’^iEj. 

3.4 Reformulation in Terms of an Eddington Ratio T 

Using the diffusion equation to relate the radiation flux Trad and the 
radiation pressure gradient, the momentum equation (eq.^ can be 
rewritten as 

V J - (£,ad - £Edd) = (T - 1) ^ (34) 

ar c r- 

where = cg/x, g is the gravitational acceleration, T = 

Frad/T^Edd. and we have neglected gas pressure. 

Equation shows that it is possible, of course, to reformu¬ 
late the results in this section in terms of an effective Eddington 
ratio T. We believe, however, that this somewhat obscures the un¬ 
derlying physics. In particular, the model developed in this section 
is most appropriate precisely when the radiation flux is small com¬ 
pared to the energy flux in the wind; and hence T are thus not 
the key dynamical variables. Moreover, we show in §4.1| that nu¬ 
merical simulations without any radiation diffusion reproduce well 
the analytic models developed here; all that is required is a fluid 
with a radiation (7 = 4/3) equation of state. And we show in ^4.2| 
that inclusion of radiation diffusion does not significantly change 
the properties of the solutions when the constraint E„ > of 
§3.3| is satisfied. The key point is that when the optical depths are 
sufficiently high, photons are trapped, the outflow is adiabatic, and 
the wind energy is roughly conserved outside the heating region. In 
our models, the outflow properties in this limit are insensitive to the 
exact opacity because once the photons are trapped out to the sonic 
point, it doesn’t matter how trapped they are. 

4 NUMERICAL WIND SOLUTIONS 

In this section we present time-dependent numerical hydrodynamic 
models that further elucidate the physics of super-Eddington winds. 
These demonstrate good agreement with the analytic solution in ^ 
We present two sets of numerical models. The first ( §4. l| l are one- 
(ID), two- (2D), and three dimensional (3D) simulations using the 
Eulerian hydrodynamics code ELASH of the model problem solved 
analytically in ^ a 7 = 4/3 fluid subjected to external energy 
input. The principal contribution of these ELASH simulations is 
that they demonstrate that multi-dimensional effects in general, and 
convection in particular, do not noticeably change the properties of 
the wind solution. The second numerical models we present are ID 
hydrodynamic simulations with MESA ( ^4.2^ . The key features of 
this simulation are that it includes radiation diffusion as well as a 
realistic stellar progenitor. 

4.1 Time-Dependent Hydrodynamics with FLASH 

In presenting our models evolved with FLASH, we focus on com¬ 
paring the global wind quantities with the analytic predictions in 
steady-state, and on the potential role of multidimensional effects. 
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Figure 4. Evolution of the spherically symmetric wind model L046R2.5-ld with FLASH, illustrating the transition to steady-state. Panel (a) shows profiles of 
density (red) and mass loss rate (blue) at the labeled times (in units of ri,lvcsc(n), rounded to two significant digits). The gray shaded region denotes the radial 
range where heating is imposed. Panel (b) shows sound speed (black) and radial velocity (green) at the labeled times. 


4.1.1 Problem Setup and Models Evolved 

We use FLASH3 jFryxell et al.|2000] |Dubey et al.|2009t to solve 
the equations of mass, momentum, and energy conservation, sub¬ 
ject to the gravity of a point mass and additional energy input. The 
public version of the code has been modified to allow for a grid 
of non-uniform spacing in spherical polar coordinates (|Femandez| 
|2012[ |Fernand^|2015^ . The equation of state is that of an ideal 
gas with adiabatic index y = 4/3. No radiation transfer is included. 
A constant heating rate per unit volume is applied in the radius 
range [r* - Ar*, r;, + Ar;,]. The FLASH models thus solve the iden¬ 
tical problem posed in ^ but include time dependence and multi¬ 
dimensional effects. 

We initialize the analytic profile in equation Q, with vanish¬ 
ing velocities. We adopt a unit system based on the radius r;,, the 
initial density p*, and escape speed v^scdh) at the location where 
heating is applied. The problem is determined by three dimension¬ 
less numbers: the initial envelope radius Rlri,, the energy deposition 
rate E Krdpkfesdfhf) and the width of the heating region 2Ar*/r;,. 
For r > R, the domain is initially filled with an ambient medium 
of constant density pa^b £ 10^'V/n so that the mass in this am¬ 
bient medium is negligible compared to that in the envelope and 
subsequent wind. The initial pressure in the ambient medium is 
Pamb = GMPmiolf- The details of the initial ambient medium at 
r>R are irrelevant to the subsequent evolution. 

The computational domain extends from r = 0.25r;, to r = 
250r/,, with logarithmic radial spacing. In 2D and 3D, the polar 
angle 6 spans the range [45°, 135°] with uniform spacing. In 3D 
the azimuthal angle (j) spans the range [0,90°], also with uniform 
spacing. The baseline resolution is such that Ar/r ^ Ad = A(l) ^ 
0.5° - 0.01 rad. The inner radial boundary condition is such that 
the ghost cells are filled with the continuation of the solution im¬ 
plied by equation 0- This forces the system to achieve steady state. 
The outer radial boundary condition is set to outflow: zero gradient 
in all variables except the velocity, which is set to zero if negative 
in the last active cell, or proportional to if positive, ensuring 
constant mass flux in the ghost cells. In 2D and 3D, the angular 
boundary conditions in {9, </>] are periodic. 

A key property of the FLASH models is that the gas in the 
polytropic atmosphere and wind is not self-gravitating. As a result, 
as noted below equation^] the Bernoulli parameter of the matter 



Figure 5. Comparison between steady-state quantities from FLASH mod¬ 
els (Table[^ gray circles, blue triangles, and black crosses) and the analytic 
wind model (red curves; eq. |24| and Fig.[^. There is essentially no differ¬ 
ence between the ID and multi-dimensional FLASH models. Top: Mass 
loss rate normalized by its critical value (eq. |23j . Middle: Wind energy loss 
rate (eq. |35^ normalized by the input heating rate. Bottom: Ratio of sonic 
point radius r, to initial envelope radius R. 
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logio radial velocity [vr/vesc(r;,)] 


-2.5 -2.0 


-1.5 


-1.0 


-0.5 



Figure 6. Radial velocity on the equatorial plane at different times in the evolution of the 3D FLASH simulation L046R2.5-3d. At early times convection 
is important at small radii, and the expansion of the wind is not completely spherical. At later times, however, when the wind has reached steady state, the 
velocity field at large radii is smooth and mostly spherically symmetric. Times are labeled in units of ri,/vescirh), rounded to two significant digits. 


interior to the heating radius is given hy Be = -GMIR. Note that 
this remains true even after the wind has reached steady state and 
the matter that was at radius r ~ R has been unbound. The reason 
is that the density, pressure, and sound speed interior to the heat¬ 
ing region do not change significantly even after the onset of the 
wind; hence Be does not change either. As a result, the critical di¬ 
mensionless number that determines the properties of the wind is 
Vesc(f^)/Vcrif This amounts to taking / = ri,/R in the analytic model 
in ^ so that f'^v^scin) = Ves^(R). 

Table [T] shows the models we evolved. We vary the heating 
rate or initial envelope radius to obtain a range of values of the 
ratio Vesc(l^)/Vcrit(£, Afenv)- Most models are run in ID and 2D to 
study the effects of convection. One model is run in 3D to quan¬ 
tify possible differences introduced by the detailed properties of the 
turbulence in the transition to steady-state. In most cases, the half¬ 
width of the heating region is Ar* = 0.075r/,; a comparison model 
with twice the width yields identical results, hence this choice is 
not important. 


4.J.2 Overview of Evolution 

Figure|^shows the evolution of the density and mass loss rate (Fig. 
1^) and velocity and sound speed (Fig.|^) in a typical spherically- 
symmetric wind model. The envelope begins to expand at radii 
r > ri, on a thermal timescale (eq. [TJ. The increase in the entropy 
caused by the injection of heat leads to a decrease in the density and 
an increase in the sound speed around the heating region. As the 
system evolves, a power-law density profile is established outside 
the heating radius. This is qualitatively similar to the MESA model 
shown in Figure but in the present case the power-law density 
profile is due to a wind rather than an extended convective enve¬ 
lope. Figure 1^ shows that the mass loss rate adjusts to a constant 
value as the envelope continues to expand, until all of the compu¬ 
tational domain reaches steady-state. The specific model in Figure 
|^(L046R2.5-ld) has a relatively high heating rate and thus at early 
times a shock forms at the outer edge of the expanding envelope. 


In the steady state wind there is no shock. In addition, for lower 
heating rates, the expansion is fully subsonic without a shock. 

The numerical solutions eventually reach a nearly perfect 
steady-state throughout the computational volume. Table shows 
the time needed to achieve steady-state in the position of the sonic 
point Tj, which we define as the first instant at which d In fj/d In t = 
10“^. This time compares favorably with the thermal time estimated 
in equation 0- For example, the thermal time and time to steady- 
state for model L005R2.5-ld are fth - 5.5 x 10^ r/,/Vesc(ni) and 
Atsteady - 1-4 X 10'' r;,/vesc(r;,), respectively. 


4.1.3 Steady-State Wind Properties &■ Comparison to Analytics 

The models shown in Table [T] cover a factor ~ 10 in the ratio 
Vesc(R)/Vcrit and straddle the point Vesc(R) = t'cnt, which marks a 
qualitative transition in the wind model of equation \2A) . We com¬ 
pare three quantities to the analytic predictions: the normalized 
mass loss rate M/Mcrit, with Merit given by equation 1 23 i, the ra¬ 
tio of wind power to input power E„IE, with 


K-f 


= I dCldpVr 


1 


r P 


j-lp 


GM 

r 


(35) 


evaluated at some r » r;,, and the ratio of the sonic point to the ini¬ 
tial envelope radius rJR, which is a proxy for via equation 1 14 1 . 


In the analytic solution, these three properties are a function of the 
ratio Vesc(R)/Vcrit only. In calculating Merit and Verit for the numerical 
models, we use the exact value of the initial envelope mass (instead 
of the approximate definition in eq.|18|l. 


M<"r’ =4;rrAV/.ln(R/r,). 


(36) 


Table [T] also reports the mass loss rate in our numerical units 
(rdphVcscdh)), which are independent of the definition of Menv 
Figure 1^ compares the steady state wind properties from the 
ID, 2D, and 3D FLASH solutions to the analytic model. There is 
good overall agreemenl in both the functional form and normal¬ 
ization of the results, suggesting that equation \2A\ is a reasonably 
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Figure 7. Comparison between the velocity fields of ID, 2D, and 3D wind 
models with FLASH. Top: Profiles of absolute value of the mean Mach 
number (eq. |39) for model L046R2.5-3d and its ID and 2D equivalent 
models (L046R2.5-ld and L046R2.5-2d, respectively). There are modest 
differences at large radii and early times, but in the steady state wind phase 
the velocity profiles are remarkably similar in all cases. Bottom: Profiles of 
mean radial velocity and r.m.s. turbulent velocity (defined as in eqns. |37|38l 
for model L046R2.5-3d at different times (as labeled in units of r/,/Vesc(r;,), 
rounded to two significant digits). The dotted lines show the turbulent ve¬ 
locity field in 2D at / = 35 and 2800. Note that the turbulent (convective) 
motions become unimportant at late times once a steady wind is established. 


good description of the system dynamics. Note that the extent of 
the quantitative agreement is somewhat sensitive to the definition 
of Menv The value adopted in equation j36| > has been used con¬ 
sistently and is the most accurate, although other definitions could 
potentially be used. 

It is evident from Figure|^and Table[^that the ID, 2D, and 
3D results are nearly identical. This suggests that multidimensional 
effects in general (and convection in particular) have no significant 
role in maintaining the wind once it has reached steady-state. This 
result is consistent with the argument given in o 

4.1.4 Multidimensional ejfects 

While multidimensional flows appear unimportant in maintaining 
the steady wind, they can make a quantitative difference in the tran¬ 
sition to steady-state. Figure [^illustrates the overall magnitude of 
multidimensional flows, showing snapshots of the radial velocity 
on the equatorial plane at a few instants in the evolution of the 
3D model L046R2.5-3d. Initially, convection is important at small 
radii, and the expansion of the wind is not completely spherical. 
At later times, however, the velocity field at large radii becomes 


increasingly smooth, with significant non-sphericities maintained 
only near the heating region. 

We quantify the effects of non-spherical flows by taking an¬ 
gular averages of velocities. The mean and r.m.s fluctuation of a 
quantity A at a given radius are 


<A> 

AA 


f dfipA 
f dQp 


(< a 2 > - (Ayf". 


(37) 

(38) 


Figure [^ shows profiles of the absolute value of the mean radial 
Mach number 


M = 


Kvr)l 

Vt<c?> 


(39) 


for model L046R2.5-2d, compared with its ID and 2D equivalents. 
At early times, the leading edge of the expanding wind lies at the 
same location in all cases, while later the 2D model evolves faster 
than 3D and ID. Nevertheless, the final Mach number profiles are 
indistinguishable once the system has reached steady-state. 

Figure [^ shows the mean radial velocity and the r.m.s fluc¬ 
tuation of the total velocity (eq.|38|l, providing further insight into 
the magnitude of multidimensional flows. At early times, convec¬ 
tive motions are strongest in regions immediately above the heat¬ 
ing radius ri„ with convective Mach numbers ~ 0.1. Convection is 
slightly more vigorous in 2D, consistent with previous numerical 
stellar convection studies (e.g., |Meakin & Arnett|2007) and as ex - 
pected from the inverse turbulent cascade in 2D ( |Kraichnan|1967| l. 

As the wind expands to larger radii, the convective motions 
subside, reaching a minimum at the heating radius, and extending 
to regions r < r;,. Once the system has reached steady-state, the 
magnitude of the turbulent velocity Av is a few % of the mean wind 
flow at the sonic point, and even weaker at larger radii. This justifies 
neglecting convection in the analytics in ^ 


4.2 Time-Dependent Hydrodynamics with MESA 

In this section we use the implicit hydrodynamics capabilities of 
the MESA stellar evolution code ( [Paxton et al.|20l51 > to calculate 
an example solution for the hydrodynamic response of a star to 
super-Eddington energy input. We utilize the analytic results from 
^and the idealized numerical experiments in §4.1| to interpret the 
more complete MESA models in this section. This calculation dif¬ 
fers from those presented in ^in that the latter were hydrostatic 
and so did not have the option of producing a wind. The MESA 
models differ from the FLASH models of the previous subsection 
in multiple ways, perhaps most importantly by including radiation 
diffusion and a realistic stellar progenitor and equation of state. The 
former allows us to explicitly check that neglecting radiation trans¬ 
port is a reasonable approximation for E„, > LEdd- 

The inlists for our MESA hydrodynamic calculation are given 
in Appendix |A2| The key properties of this model include that 
we utilize MESA’s implicit hydrodynamics capabilities, which 
amounts to solving the time dependent spherically symmetric hy¬ 
drodynamics equations. This includes radiation in the diffusion ap¬ 
proximation assuming local thermal equilibrium. We do not utilize 
the MLT-l-l- convection module used in ^ since this would over¬ 
estimate the efficacy of convective energy transport. We also turn 
off MESA’s wind models so that the mass on the domain is con¬ 
served. The one significant simplification in our calculation is that 
we set the opacity to be electron scattering everywhere. Calcula¬ 
tions with a constant opacity 100 times larger yielded very similar 
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Figure 8. Density and mass loss rate M = Anr^pv profiles in hydrodynamic MESA calculations, at dift'erent times after the onset of energy deposition with 
E = 10^Lq at r/, ^ 35Rq (in a 23 Mq, R = 60Rq, Z = Zq model at core He exhaustion). The initial thermal time (eq.^ at /*/, is 0.4 yr, which sets the initial 
expansion time of the envelope. At late times the heating has generated a quasi-steady state wind with M ^ 0.4M© yr“^ (independent of radius over a factor 
of ~ 100 in radius). In addition, the higher speed wind at late times shocks on the slow wind at earlier times (see the velocity profiles in Fig.[9j, generating a 
swept up shell of mass visible as the excess density/M at large radii. 


results for the wind properties, demonstrating that radiation diffu¬ 
sion is not critical for the outflow properties, consistent with the 
arguments in 0 Our restriction to a constant opacity is primarily 
because of complications associated with the precipitous drop in 
opacity when hydrogen recombines below < 10'' K. This recom¬ 
bination is likely important for understanding the outer portions of 
super-Eddington stellar winds, including the observational signa¬ 
tures of such winds. We leave a study of this important aspect of 
the problem to future work (but see ^for a brief discussion). 

We carried out calculations with a variety of massive stel¬ 
lar progenitors and energy deposition rates/locations. We present 
one illustrative example here. We consider a 30 Mq, Z = Zq 
stellar model evolved to He exhaustion, at which point its mass 
and radius were 23Mq and 60Rq, respectively|^We then deposited 
E = 10^ Lq in a region centered at rn = 35Rq. Given the density 
profile prior to heating in this stellar model, the dimensionless en¬ 
ergy injection rate as defined for the FLASH models in §4.1| is 
E ! diPdhtvescdhf - 6 X 10“^. This is similar to the heating rate 
in the FLASH model L046R2.5 (ID, 2D, and 3D) shown in Fig¬ 
ures EE As time goes on, however, p(r;,) decreases in the 

MESA models (as discussed below) so that the dimensionless E in¬ 
creases, making the MESA model at late times more analogous to 
the higher E FLASH models. 

The left panel of Figurej^shows the density profile and mass 
loss rate M = 4ndpv at early times in the response to energy depo¬ 
sition in the MESA models while the right panel of Figurej^shows 

^ These stellar parameters differ from the models in despite being for 
the same ZAMS stellar properties and stellar wind model. This is in part 
because here we use electron scattering opacity only. Moreover, not us¬ 
ing MLT-I--I- changes the stellar structure and the integrated effect of stellar 
winds on the model by the time of He exhaustion. 


the same quantities during the phase when a strong steady wind is 
established. Figure|^shows the velocity profile at these same times. 

The total stellar mass initially exterior to the heating radius is 
~ 0.3Mq so that the initial thermal time of the envelope is ~ 0.4 
yr. The stellar envelope initially expands outwards on of order the 
initial thermal time, driven by the excess thermal energy deposited 
near ~ r*. The dynamics in this phase is reminiscent of the hydro¬ 
static models in ^but even at this early stage the density profile is 
somewhat flatter than the convective models in Figure [T] The evo¬ 
lution of the density profile in the MESA model is very similar to 
that seen in the FLASH simulations in Figure]^ The mass outflow 
rate in Figure[^is relatively independent of time at == O.4M0yr“*, 
even at early times, but the velocity of the flow continues to accel¬ 
erate reaching ~ 300 km s“' at f = 6 yr. A quasi-steady outflow 
only develops around ? = 3 yr. As a result, at early times M is best 
interpreted as being due to the nearly hydrostatic expansion of the 
stellar envelope, rather than a wind. During the wind phase at t > 4 
yr, M becomes nearly independent of radius over a wide range of 
radii, as expected for a steady wind. 

The increase in the wind velocity with time seen in Figure 
1^ leads to a pile up of mass at large radii as the wind runs into 
mass ejected at earlier times. This accounts for the shell of matter 
at radii ~ lO^-^-^'Ro (Fig. § and the strong shock evident in the 
velocity profile at late times (Fig.|^. At this stage of evolution, the 
strong stellar wind is effectively inflating a ‘wind bubble’ inside the 
remnants of the outer stellar envelope. This is not as visible in the 
FLASH simulations in §4. l| in part because FLASH is an Eulerian 
code and thus matter at large radii leaves the computational domain. 

Figure [T^ shows the energetics of the outflow at f = 6 yr. 
The ‘wind’ luminosity is defined as E„ = Be M while the ‘total’ 
luminosity is the sum of the wind, photon, and convective contri- 
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Figure 9. Outflow velocity and sound speed as a function of radius in hy¬ 
drodynamic MESA calculations, at different times after the onset of energy 
deposition with E = 10^L© at r/, ^ 357?© (in a 23 M©, R = 607?©, Z = Zq 
model at core He exhaustion). As the energy deposition ejects the stellar en¬ 
velope to large radii, the velocity of the outflow accelerates to ~ 300 km s“^ 
and the sonic point moves in to small radii. At late times most of the energy 
deposited in the stellar envelope is carried to large radii by the wind (see 
Fig.[^. 


butions|^The total outflow power is nearly independent of radius, 
as expected for a quasi-steady wind. Moreover, most of the energy 
flux is carried by the wind, with photon diffusion contributing a 
total luminosity =: 1.5 10 ®Lq and convection much less. The sub¬ 
dominance of convection in Figure [T^ is consistent with the ana¬ 
lytic arguments in 5 3.1 The photon luminosity in the late steady 
wind stages exceeds the Eddington luminosity by a factor of == 2, as 
is perfectly possible for an outflow (see At earlier times when 
the wind is much lower velocity, the photon luminosity is almost 
exactly the Eddington luminosity, confirming that the moderately 
super-Eddington photon luminosity at late times is due to the wind. 
Finally, note that the total wind power is about a factor of ~ 2 less 
than the energy input rate E = 10^L©. Some of the energy is lost to 
work done against gravity escaping the stellar potential. 

The analytic models in ^can explain many of the broad prop¬ 
erties of the MESA simulation. For concreteness, consider f ~ 5 yr 
when Vcrit - 140km s“^ - 470 km s“\ and M/A/env - 200. 

The numerical solution at this time has Vco - 160 km s“^ and M ^ 
0.4 M© yr“^. The model in ^predicts M ^ O.OSMci-it - 0.5 M© yr“^ 
and Voo - Vcrit - 140 km s , which are comparable to the numer¬ 
ical values. Physically, the mass loss rate corresponds roughly to 


^ To calculate Be in MESA, we calculate the gravitational potential by in¬ 
tegrating V0 = -GMrIr- in run_star.extras and add this to the enthalpy per 
unit mass (the sum of the internal energy per unit mass and P/p) and the 
kineti c ene rgy per unit mass. 

^ In j3.l| we argued that the outflow should not be convective because it 
is adiabatic outside the heating region. In the MESA calculations, photon 
diffusion means that the outflow is not strictly adiabatic, which leads to 
modest energy transport by convection. 


Figure 10. Total, wind, photon, and convective luminosities as a function of 
radius in hydrodynamic MESA calculations, 6 yr after the onset of energy 
deposition with E = 10^E© at ri, ^ 357?© (in a 23 M©, R = 607?©, Z — Zq 
model at core He exhaustion). The wind power is E,^, - BeM, where Be is 
the Bernoulli parameter. The total power is the sum of the other three con¬ 
tributions. The wind is the dominant energy transport mechanism to large 
radii at these late times when the outflow speed is ~ 300 km and the 
sonic point of the wind has moved in to small radii (see Fig.[^. Convec¬ 
tive energy transport is negligible, in contrast to the hydrostatic models in 
Figure]^ Note also that the emergent photon luminosity at the outer edge 
of the wind is larger than the photon luminosity in the wind at intermediate 
radii. This is due to the strong shock produced as the high speed wind at 
small radii encounters the slower moving shell at large radii, thermalizing 
the kinetic energy of the wind (Fig. I3&^. 


E 0.5vi:sc(f*)^V (Regime 2 in Note that in this regime we 
expect the mass loss rate M to be relatively independent of the 
asymptotic velocity Voa- This is why the mass loss rate does not 
change significantly in time (Fig.[^ even as the asymptotic veloc¬ 
ity increases and the sonic point moves in from large radii (Fig.|^. 
Figure (left panel) shows that this characteristic mass-loss rate 
holds even when there is no super-sonic outflow. This is because 
even at the early times shown here, neither photons nor convection 
carry a significant fraction of E outwards, so the majority of the 
energy input goes into lifting matter out to large radii. 

One aspect of the MESA simulations not captured by the 
steady state wind models in ^and §4.1| is the continued evolu¬ 
tion of the wind to higher speed and thus larger wind power even 
when the sonic point has moved well interior to the outer radius of 
the ejected mass (and so the outflow is now out of causal contact 
with the previously ejected matter). There is never a true steady 
state. This is not due to radiation diffusion modifying the structure 
of the envelope since the calculations with 100 times larger opac¬ 
ity yielded the same behavior. We interpret this as due to changes 
in the stellar mass and envelope structure due to the outflow (the 
total mass ejected during the simulation is =: 2.1 of order 10% 
of the initial mass). The decreasing stellar mass and the ejection of 
the stellar envelope decrease the scale height of the envelope at the 
fixed heating radius r;,. This in turn decreases Me„v and increases 
Vcrit, leading to a smaller sonic point radius and a more powerful 
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outflow at late times. It is worth stressing, however, that in reality 
the longer timescale evolution of the wind properties will depend 
in detail on how the source of heating adjusts as the envelope is 
ejected (there is no such back reaction in our MESA calculations). 


5 OBSERVATIONAL SIGNATURES 

In this section we present estimates of the observational character¬ 
istics of the wind models described earlier in this paper. For con¬ 
creteness, we assume in this section that the wind is in Regime 1 
of ^ and use the analytic estimates of the wind properties to es¬ 
timate its radiative properties. Recall that in this regime £■„, ~ E, 
i.e., the kinetic power in the wind is comparable to the total power 
E supplied to the stellar envelope (eq. [28| >. The outflow in gen¬ 
eral consists of two components ( §4.2| (: the steady wind at small 
radii surrounded by a shell at larger radii comprised of previously 
ejected matter swept up by the wind. We first estimate the observa¬ 
tional properties of the steady wind neglecting the outer shell and 
then briefly discuss the impact of the outer shell. 


5.1 The Eree Wind 


As the wind accelerates outwards at the sonic point ~ r^, photon 
diffusion is initially negligible and the trapped photon energy is 
instead carried outwards by advection in the wind ( §3.3| >. However, 
the ratio of the diffusion time dpKjc to the advection time 

4dv - r/v is oc 1/r for constant opacity in a steady wind. Thus 
at the ‘diffusion radius’ rj, photon diffusion becomes important. 
The diffusion radius can be estimated by using the fact that at r„ 
tdiff/fadv ~ El Lem, which yields 


rd ~ 



(40) 


The total bolometric luminosity radiated by the wind is set by the 
thermal energy content at and is given by 

/ E 

iwind ~ f'Edd 7 - ■ (41) 

\ r^Edd / 

Equation shows that radiated power of the wind is super- 
Eddington and increases with increasing heating rate oc Pq^ 
the MESA model in ^4.2[ equation |4 1 1 predicts L„ind ~ S.VLedd, 
which is similar to the factor of 2 found in the numerical solution. 
Note that if density inhomogeneities in the stellar envelope and/or 
outflow increase the rate of photon diffusion (e.g., |Shaviv ||200T1 
lOwocki ef al.|2004) , fhis effectively increases LEdd in equationpl] 
thus increasing the total power radiated by the outflow. 

For Ej < r < Erf, the outflow is adiabatic and so the temperature 
is given by T oc oc As a result, the temperature at is 
given by T(rd) ^ T{rP{rJrdf^^, i.e., 

T(rd) ~ 3 X 10^ K (42) 

where a:o .4 = a:/ 0.4 cm^ g“' is the opacity at Vd scaled to the electron 
scattering opacity. The temperature in equation |42| depends only 
weakly on model parameters and is sufficiently hot that the elec¬ 
tron scattering opacity is likely a reasonable first approximation for 
estimating and L„ind in equations and For more quanti¬ 
tatively accurate estimates the diffusion radius can be obfained 
from tdiff - 4dv using a more realistic opacity K{p, T). 

The effective temperature associated with the observable 
emission L„ind will in general be less than T{rd) in equation[4^be- 
cause the optical depth at Vd is still quite large ~ c/v„. Formally in 


the free steady wind, the photosphere is at a very large radius such 
that the time for the outflow to reach the photosphere is hundreds of 
years for our fiducial massive star parameters. This implies that in 
practical cases of interest, the photosphere is limited by the radius 
the wind has reached at time t, r ~ Vct. This yields an effective 
temperature of 

( ,, v-l/lO 

lo^j 


where tyr is the time since the steady outflow was initiated in years. 

The applicability of equation|^is limited by several factors. 
Once Ted falls below ~ 5000 K (i.e., after a few months), recombi¬ 
nation is important and will fix the effective temperature to ~ 5000 
K. The photosphere will then lie inside the outer radius of the wind 
that is continuing to move outwards. An additional complication is 
that stars with mass loss rates a s large as equatio n |25| can form dust 
outside a radius of ~ 10*^ cm i Kochanek|201 l| l. This implies that 
after of order a few -10 years dust will begin to form and will repro¬ 
cess some of the wind emission to even longer wavelengths. Prior 
to efficient dust formation, however, the estimates of this section 
demonstrate that the super-Eddington wind is likely to be a bright 
super-Eddington optical source. At later times the emission will be 
increasingly in the infrared. 


5.2 The Swept-Up Shell 

The estimates of the previous section assume that once the wind is 
accelerated out through the sonic point that it continues to expand 
unimpeded to larger radii. This is in fact not true in the MESA 
models in §4.2| (Figs and [^. Instead, at early times most of the 
energy input goes into slowly expanding the atmosphere outwards. 
At later times, the density in the heating region is lower so the wind 
can be accelerated to higher speeds. The end result is a compara¬ 
tively slowly moving shell of matter at larger radii which confines 
(in ID) the higher speed wind at late times. The high speed wind 
shocks on the outer shell, converting its kinetic energy into thermal 
energy. 

For a shell of mass at radius r, the optical depth through 
the shell is ~ 60(Afsheii/MQ)(r/10*^ cm)“^^(:o 4 . At large radii dust 
will form (e.g., |Kochanekp011^ at which point the optical depth 
will be even larger than this estimate. Thus much of the shocked 
wind energy may be thermalized and re-radiated ( |Smith|2013| ar- 
gues for a model along these lines for Eta Carinae, in which much 
of the emission of the great eruption is powered by circumstellar 
interaction). This implies that in many cases the radiated power by 
super-Eddington winds is likely to exceed the estimate in equation 
l^and be closer to the kinetic power in the wind. Note that this is 
true in the MESA model shown in Figure[T^ where the emergent 
photon luminosity at the outer edge of the wind is larger than the 
photon luminosity advected out in the bulk of the wind by a factor 
of - 2. In detail, the efficiency of this “internal shock” emission 
will depend on how the geometry and kinematics of the wind and 
the swept-up shell at large radii evolve in time. 


6 APPLICATION TO MASSIVE STARS 

We now briefly discuss the application of our results in the previous 
sections to outflows from massive stars. To provide some context 
for this application. Figure [TT] (left panel) shows the thermal time 
(defined as in eq.^ as a function of envelope mass for five different 
massive star models, from compact Wolf-Rayet-like models (R — 
0.7Rq) to red super giants (R == 10 ^ Rq). The models in Figure[TT]are 
all at the end of He core fusion (see Appendix|A3[). Recall that the 
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Figure 11. Left: Thermal time for the stellar envelope (eq.^ as a function of envelope mass M(> f) (mass exterior to a given radius r) for different massive 
Stellar progenitors at He core exhaustion: (M/Mo,R/^o) = (9.3,0.76), (11.1,2.6), (11.3,18.7), (11.43,275), and (19.4, 1100), from left to right in the plot, 
respectively. The stellar models are labeled roughly by whether they would be classified as a Wolf-Rayet star (WR), a blue supergiant (BSG), or a red 
supergiant (RSG). The thermal time sets the timescale for the envelope to respond to energy injection, and thus the time for the wind solutions calculated in 
this paper to apply. It is scaled here to a heating rate of 10^L©, which is roughly 30 times the electron-scattering Eddington luminosity for a 10 Mq star. Mass 
ejection via super-Eddington winds on decade (e.g, Eta Carinae) or month-decade (e.g., Type Iln) timescales requires ene rgy deposition in the outer ~ 0.1 - 1 % 
of the star by mass, with the exception of the largest RSGs. Right: Dimensionless parameter / = 2Be(r)}vl^^{r) (eq. |l0) that characterizes the binding energy 
per unit mass of the stellar envelope, for the same models as in the left panel. Lower values of / in general lead to more powerful winds (^& Fig.[^. Stellar 
models with larger envelopes (BSG and RSG) have lower values of / and are thus more prone to high speed energetic winds. 


thermal time scale at the heating radius sets the overall timescale 
on which the stellar envelope adjusts to the additional energy input 
and thus the timescale on which the envelope expands and a steady 
wind develops. Note that for the compact progenitors in Figure [TT] 
hhemiai is simply oc M(> r) because the radius r does not change 
significantly for the envelope masses considered. This is not true 
for the extended progenitors. 

The right panel of Figure also shows the dimensionless 
binding energy / = 2Be{r)lvl^^{r) (eq. 10 1 of the stellar envelope. 


For everything else fixed, lower values of / lead to more energetic 
winds (^& Fig.[^. Figure [TT] shows that the stellar models with 
extended envelopes (BSG and RSG) have lower values of / and are 
thus more prone to high speed energetic winds. This is qualitatively 
consistent with the polytropic atmosphere calculation in §3.2| 


6.1 Eta Carinae and Luminous Blue Variables 

LB Vs are the most dramatic manifestation of episodic mass loss in 
massive stars (see, e.g., [Davidson & Humphreys|2012[[Smith|2014| 
for reviews). The prototypical (albeit extreme) example is Eta Cari¬ 
nae whose great eruption from ~ 1840 - 1850 lasted over a decade 
with the photon luminosity exceeding ~ 10 times the Eddington 
limit for a 50 Mq star throughout this period. The properties of the 
gas ejected during the great eruption can be estimated from the sur¬ 
rounding nebula, which yields a mass, velocity, and kinetic energy 


of ~ 20Mq, ~ 500km s , and ~ 5 lO'^^erg, respectively (Smith 


|et al.||200^ . The total photon energy radiated during this period 
was of order, though probably somewhat less than, the kinetic en 
ergy of the ejecta. 


The physical mechanism responsible for LBV outbursts re¬ 
mains uncertain. Nonetheless, the existence of super-Eddington ra¬ 
diative and kinetic luminosities for many dynamical times suggests 
that the wind models developed in this paper should be reason¬ 
ably applicable. Our models clarify the conditions required for the 
wind kinetic energy flux to be of order the energy input rate E 
and larger than the radiated photon luminosity. This condition must 
roughly be satisfied in LB Vs given that the energy requirements to 
explain the ~ 10^° ergs of kinetic energy in the outflow would go 
up significantly if E„ « E. Our results show that Ew E > Trad 
if the critical speed Verjt (eq. |22| l is of order the escape speed at 
the heating location. When this condition is satisfied, the outflow 
velocity is naturally a few Vcru ~ 300km s“' (eq. |^& Figure]^ 
and can be as large as that of the escape speed from the heating 
region. If the current stellar mass and radius of Eta Carinae reflect 
the structure prior to the great eruption, Vesc(r*) ~ 400km s“* if 
the heating occurs near the surface, as is required to have a short 
thermal time (Fig. [Uj. Thus, the characteristic velocities predicted 
by our model are of order that observed for much of the mass in 
the great eruption. We suggest that the most likely explanation for 
Vcrit ~ Vescfr*), and hence an efficient conversion of energy input to 
wind kinetic energy, is that the overall evolution of the density pro¬ 
file in response to mass loss is qualitatively similar to that found in 
our MESA calculations in ^4.2| (Fig.[8t. In particular, if the typical 
mass in the stellar envelope decreases in response to mass loss, the 
star will naturally approach the limit in which much of the mass is 
ejected with E„, ~ E. 

The month-decade timescale of LBV eruptions favors signif- 
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leant energy deposition relatively close to the stellar photosphere. 
In particular, the thermal time of the stellar envelope is likely < a 
few years to explain this evolutionary timescale. For most stellar 
progenitors this requires that the excess energy input occur in the 
outer few percent of the stellar mass (Fig. [nj. It is important to 
stress that this timescale constraint on the depth of the heating re¬ 
gion does not limit the mass ejected to be less than a few percent of 
the stellar mass, because as the outer stellar envelope is shed, new 
matter from deeper down in the star replenishes the outer envelope. 


6.2 Pre-Supernova Mass Loss & Type Iln Supernovae 


Up to ~ 10% of supernova (SN) progenitors appear to experience 
episodes of intense mass loss in the century to weeks leading up to 
core collapse (e.g., [Smith et al.|201 l[|Ko c hanek|201 l[|Kiewe et al.| 
|2012[ see |Smith|2014 for a review). We have previously proposed 
that powerful mass loss in the last year-decade prior to core col¬ 
lapse can be produced by vigorous convection in the stellar core 
exciting a super-Eddington wave flux that travels out into the stel¬ 
lar envelope and dissipates jQuataert & Shiode||2012| jShiode &| 
|Quataert||2014| >. In addition to directly driving mass loss, wave- 
driven inflation of the stellar envelope could also trigger mass loss 
via Roche Lobe overflow in close binary systems JSmith & Amettj 
[M^pcley & SokerpOll} . 


The wind models developed in this paper are particularly ap¬ 
plicable to the wave-driven mass loss mechanism since it directly 
corresponds to a non-local redistribution of energy in the star that 
can lead to super-Eddington heating rates in the stellar envelope. 
Our wind models can explain the characteristic mass loss rates of 
~ 0.01 - 1 Moyr^' and velocities of ~ 200 - 500 km s“* inferred 
from many Type Iln SNe (e.g., [Kiewe et al.|2012} , if a reasonable 
fraction of the wind is ejected roughly in ‘Regime F of ^ in which 
the wind energy flux is of order the energy input rate E (which is 
the wave energy flux from the stellar core in the present context)]^ 
One seemingly fine tuned aspect of this model is the need for 
energy to be deposited close the stellar surface in order for the ther¬ 
mal time to be short enough (Fig. [nj to explain enhanced mass 
loss in the months-year prior to core collapse (the same concern of 
course arises in other models as well). In fact, however, this is a nat¬ 
ural property of the wave-driven mass loss model, at least for some 
progenitors. The reason is that the outgoing waves in the stellar en¬ 
velope are sound waves that damp primarily by radiative diffusion. 
Thus they necessarily only dissipate their energy relatively close 
to the surface where the thermal time across the wavelength of the 
waves is shorter than the group travel time. 


7 DISCUSSION 

In this paper we have studied the physical properties of super- 
Eddington stellar winds. By this we specifically mean winds in 
which the kinetic power approaches or exceeds the Eddington lu¬ 
minosity. This work is motivated by phases in stellar evolution in 
which super-Eddington energy deposition can heat a region near the 
stellar surface, potentially driving a powerful wind. Examples of 
this phase likely include classical novae, where the resulting mass 


^ In some Type Iln SNe, the wind speed inferred from the narrow hydrogen 
lines approaches 1000 km s“*. In our model this requires compact progen¬ 
itors with a large escape speed and wave heating at relatively low densities, 
i.e., close to the photosphere. E.g., we find wind speeds of this magnitude 
in MESA calculations with compact BSG progenitors. 


loss can limit the ability of white dwarfs to approach the Chan¬ 
drasekhar mass; radius expansion X-ray bursts, whose properties 
are i mportant for constraining the neutron star equation of state 
(e.g., Ozel||20d^ ; luminous blue variables (LBVs), whose year- 
decade long outbursts may dominate the mass loss from massive 
stars (e.g., [Smith & Owocki||2006^ ; and enhanced pre-supemova 
mass loss inferred via circumstellar interaction in Type Iln super¬ 
novae. The sources powering this excess energy deposition are var¬ 
ied, potentially including unstable fusion, wave redistribution of 
energy, large changes in opacity with temperature, or the effects of 
a companion in a binary system (e.g.,[Schwarzschild^Harm|jj^ 


[Paczynski|1976[[Quataert & Shiode|2012[ Smith & Amett|2014^ . 

When the energy flux in a stellar wind is super-Eddington, 
the photons in the resulting outflow are trapped with a diffusion 
time long compared to the outwards advection time ( §3.3[ l. In this 
limit the outflow is essentially adiabatic in the key region where 
the wind mass loss rate, velocity, etc. are set (exterior to the heat¬ 
ing region, but interior to the sonic point). Moreover, we have 
shown analytically and numerically that convection is unimportant 
in quasi-steady super-Eddington stellar winds ( |3.1[ & [4.1[ l. This 
is in contrast to the roughly hydrostatic response of stars to super- 
Eddington heating, which inevitably drives convection that rear¬ 
ranges the structure of the star (e.g., [Josset al.|197^ see 

The fact that photons are trapped and convection is sub¬ 
dominant simplifies the dynamics of super-Eddington stellar winds. 
It implies that they can be modeled as a y = 4/3 radiation domi¬ 
nated fluid together with a physical prescription for the source of 
additional heating. This is a version of [Parker] ( [1958[ l’s solar wind 
model generalized to a radiation dominated fluid. The fact that the 
photons are trapped also argues against radiation diffusion through 
a porous stellar atmosphere as being the key physics that sets the 
wind properties (as was suggested in previous work; e.g., |Owocki [ 
[et al.|2004l l: the rate of photon diffusion must be greatly enhanced 
for the wind to become non-adiabatic and thus for the wind prop¬ 
erties to be significantly modified from those calculated here. 

We have analytically solved for the properties of super- 
Eddington winds for a simple model in which heating at a rate E oc¬ 
curs primarily near a heating radius r;,. The wind properties depend 
on the ratio of the escape speed at the heating radius Vesc(r*) to a 
characteristic speed in the problem Vent = (see Fig. 

where is the fraction of the stellar mass near the heating 

radius. For Veru ^ Vesefr*), the wind kinetic power is of order the en¬ 
ergy input rate E. The asymptotic wind speed is ~ 300 km s“* for 
typical parameters relevant to massive stars (eq.[26[> and is bounded 
from above by Vesc(r;,). In this regime, the outflow is related to what 
is sometimes termed the “photon tired limit,” typically taken to be 
when most of the photon luminosity of the star is converted into 
wind kinetic energy (e.g., [Owocki & Gayle^[1997[ ). In fact, the 
models developed here can exceed this nominal limit on the wind 
kinetic power. The reason is that we do not assume that the heating 
rate in the stellar envelope is bounded by the photon luminosity. 
This is indeed not necessarily the case for unstable thermonuclear 
fusion, waves excited in the stellar core heating the envelope, or 
heating due to an external companion (e.g., tides, common enve¬ 
lope). In our models, the maximum wind power is instead bounded 
by the energy input rate E into the stellar envelope. 

For Verit £ Vesc(r;,), the asymptotic kinetic power in the wind 
is less than the excess energy supplied to the stellar envelope E. 
This is because most of tbe energy goes into work against gravity 
unbinding mass from the stellar potential. Indeed, the stellar mass 
loss rate in this regime is set by £ == 0.5MVesc(r;,)^ (eq. |30f . The 
asymptotic speed is typically < Verit (eq.j^and Fig.|^. 
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We have validated the analytic model developed here using 
hydrodynamic simulations with FLASH and MESA (the lat¬ 
ter utilize the new hydrodynamic capabilities of MESA; |Paxton| 
|et al.|20f5) . The purpose of including simulations with two differ¬ 
ent codes is that each explicitly tests different assumptions of the 
analytic model. In particular, the FLASH simulations use a simple 
7 = 4/3 equation of state and do not include radiation diffusion 
but are multi-dimensional while the MESA simulations are one¬ 
dimensional but include radiation diffusion and a realisfic equation 
of state and stellar progenitor. 

The key result of the ELASH simulations is that the steady 
state solutions of the model analytic problem in multiple dimen¬ 
sions (largely 2D, but one 3D simulation) agree remarkably well 
with the analytic results, and with one-dimensional simulations 
(Fig. s . This is fundamentally because convection can be important 
in the initial hydrostatic expansion of a stellar envelope in response 
to heating, but it is unimportant in the resulting steady state wind 
that develops (see Fig § 

Our MESA hydrodynamic simulations model super- 
Eddington energy input into the envelope of a massive star. During 
the first thermal time, the envelope mass at the heating location is 
so large that there is not a strong supersonic wind: the envelope is 
essentially hydrostatically lifted to large radii (Eigs.[^&|^. As the 
mass in the stellar envelope decreases, however, a supersonic wind 
quickly forms and accelerates outwards, effectively inflating a 
stellar wind bubble inside the slower matter ejected at earlier times. 
The wind in this phase has a mass loss rate and velocity similar 
to that predicted by the analytic models in 0 This explicitly 
demonstrates that radiation diffusion does not strongly affect the 
wind properties when the kinetic power in the wind is of order or 
larger than the Eddington luminosity. 

Although the calculations we have presented are general and 
potentially applicable to a range of astrophysical environments, we 
have in particular highlighted their application to powerful outflows 
from massive stars (^. This includes both LBVs outbursts such 
as Eta Carinae’s great eruption and the large mass loss rates that 
precede some core-collapse supemovae. 


LBV outbursts are characterized by a super-Eddington photon 
luminosity and an outflow with a super-Eddington kinetic power 
(e.g., [Smith et al.|2003| l. The fact that this activity is on for many 
dynamical times argues against a pure shock-mediated phenomena 
and in favor of an energy source that drives a continuous wind. 
Although the ultimate energy source for LBV eruptions is poorly 
understood, we argue that our models can explain the characteris¬ 
tic mass-loss rates of ~ 1 Mq yr“’ and ve locit ies of several hun¬ 
dred km s“* of much of the ejected mass ( ^6.l| l. We speculate that 
the evolution during an eruption shares some qualitative similarities 
with the MESA models in Eigures [^ & |^ at early times much of 
the input energy goes into inflating the stellar envelope. Only then 
can the outflow be accelerated to high speeds. 


Our super-Eddington wind models produce a super-Eddington 
photon luminosity like that observed in LBV eruptions and some 
pre-SN outbursts from massive stars (e.g., [Smith et al.|201 l[|Of^ 
[et al.|20l^ . The luminosity of the steady wind freely expanding 
into a vacuum is given by ~ LEdd(E/LEdd)F^ (eq. 41 1 . However, 
in our MESA calculations, the strong shock driven as the higher 
speed wind at late times encounters the slower, previously ejected 
envelope can thermalize much of the wind kinetic power. Depend¬ 
ing on the exact wind kinematics as a function of time, this may 
substantially increase the photon luminosity of the wind ( ^5.2[ & 
Eig.|10[l. More realistic radiation transfer calculations of this pro¬ 


cess would be particularly valuable in quantitatively connecting the 
models developed here to observations. 

Throughout this work we have emphasized the essentially hy¬ 
drodynamic character of super-Eddington stellar winds, in which 
radiation transfer is not dynamically that important because the 
photons are trapped and advected with the fluid. We acknowledge 
that although this is true in the models we have developed, it is 
possible that more realistic calculations will show that the out¬ 
flow is sufficiently inhomogeneous to substantially increase the 
rate at which photons diffuse through the outflow (e.g., [Shaviv[ 
|2001[[Owocki et al.|2004[ >. Understanding this more quantitatively 
will ultimately require multi-dimensional radiation (magneto)- 
hydrodynamical simulations (e.g., [Jiang et al.[2015[ ). 

One particularly interesting extension of the work described 
in this paper is to the possibility of outflows generated when large 
changes in opacity with temperature near the surface of a star cause 
the luminosity to suddenly be super-Eddington. If this occurs suf¬ 
ficiently close to the photosphere, convection may be unable to 
carry outwards the energy generated in the stellar interior. In hy¬ 
drostatic stellar models, this leads to gas pressure and density inver¬ 
sions (e.g., [Joss et al.|1973[ [Paxton et al.|2013[. In hydrodynamic 
models, however, a wind may develop (e.g., Kato & Hachisu |1994[ 
[Eichler et al.|1995) . Eor this application, our calculations need to 
be extended to take into account the temperature dependence of the 
opacity and the fact that the flux may be super-Eddington over only 
a modest range in temperature (e.g., at the iron opacity bump). 

In future work it would also be valuable to extend our analysis 
to lower mass outflow rates. In particular, when the kinetic power 
of the wind is below the Eddington luminosity, it is no longer cor¬ 
rect to assume - as we have done - that the wind is adiabatic be¬ 
tween the heating region and the sonic point ( §3.3[ l. Understanding 
outflows with lower mass loss rates thus requires a more careful 
treatment of the effects of photon diffusion. 
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APPENDIX A: STELLAR MODELS 

Our stellar models are constructed using the MESA stellar evolu¬ 
tion code, version 7664. This version includes several additional 
boundary conditions for the full hydrodynamic evolution that were 
useful for the wind solution in ^4.2| (these boundary conditions are 
not in MESA release 7624). The models use the following inlist 
controls file. All runs use inlist_massive_defaults in addition to the 
specific flags given below. 

A1 Hydrostatic Models With Heating Erom ^ 

We first run a model to He core exhaustion using 
&star_j ob 

ere ate _pre .main .sequence .model = .true, 
(fecontrols 

okay .t o .re duce .gr adT .exces s = .true. 

initial.mass = 30 
Zbase = 0.02 

varcontrol.target = ld-3 

Dutch_scaling_factor = 1.9 

Dutch.wind.lowT.scheme = ‘ de Jager ’ 

Hot.wind.scheme = ‘Dutch’ 

Xa_centra 1.1 ower.1 imit_species (1) = ‘he4’ 

Xa.central.1 ower. 1 imit (1) = ld-5 

The wind parameter Dutch.scalingYactor was chosen as above to 
end up with a compact stellar model at He core exhaustion, since 
this most effectively demonstrates the effects of extra energy de¬ 
position. We then used the output model above as an input model 
and evolved it using the extra.energy subroutine in run_star_extras.f 
(described below in §A4| l with the same &controls as above except 

(festar.j ob 

set.initial.age = .true, 
initial.age = 0 

set.initial.dt = .true, 
years.for.initial.dt = 0.001 

change.V.flag = . true . 
c h a n g e . i n i t i a 1. V .f 1 a g = .true, 
new.V.flag = . true . 

(fecontrols 

use.other.energy = .true, 
x.ctrl (l) = 3d6 
x.ctrl (2) = 2. 
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x.ctrl (3) = 0.3 
X -Ctrl (4) = 0.1 
X -Ctrl (5)= 1 0. 

Varco ntro 1-target = ld-4 

Dutch-scaling-factor = 0. 

inax_age = 100. 

The choice of max_age is arbitrary. The models were stopped 
roughly after the times shown in Figures|T|&[^ 

A2 Hydrodynamic Models with Heating From ^4.2| 

We first run a model to H core exhaustion using the same inlist 
parameters as in the first model of Appendix |A1 [ except with 

&controls 

okay _to_reduce_gradT_excess = .false. 

mixing_length_alpha = 1.89 

The choice of a different mixing Jength.alpha is unimportant. We 
then run a model to He core exhaustion using 

&star_j ob 

change.v.flag = . true . 
change_initial_v_flag = .true. 
new-V-flag = . true . 

&controls 

okay _to_reduce_gradT_excess = .false. 

initial-mass = 30 
Zbase = 0.02 

V ar c ontr o 1 _t ar g e t = ld-3 

mixing_length_alpha = 1.89 

Dutch-scaling-factor = 1.9 
Dutch-wind-lowT-Scheme = ‘ de Jager ’ 
Hot-wind-Scheme = ‘Dutch’ 

Xa_central-1 ower-1 imit_species (1) = ‘he4’ 

Xa_central-1 ower_ 1 imit (1) = ld-5 

use_simple_es_for_kap = .true. 

The key differences relative to the model in Appendix |A1| are the 
use of electron scattering opacity only and the absence of MLT++. 
Use of MLT++ in the hydrodynamics calculation would artificially 
suppress the generation of a wind. For consistency, we thus evolved 
the progenitor model prior to energy deposition without MLT++. 
We then used the model above as an input model and excised the 
core using 

&star_j ob 

remove_initial_center_by_mass_Msun = 13 
&controls 

max_model_number = 1 

The removal of the core allows the subsequent hydrodynamic cal¬ 
culation to focus on the envelope where the energy deposition oc¬ 
curs. Finally, we used the output model above with an excised 


core as an input model and evolved it using MESA’s implicit 
hydrodynamics solver. We used the extra.energy subroutine in 
run_star-extras.f (described below in §A4^ with the same &controls 
as above except 

&star_j ob 

relax_to_this_tau_factor = ld-4 
d1ogtau_factor = .1 
r e 1 a X _ i n i t i a 1 _ t a u _ f a c t o r = .true. 

change.E.flag = .false. 
new-E-flag = .false. 

s e t _ i n i t i a 1 _ a g e = .true, 
initial.age = 0 

s e t _ i n i t i a 1 _ d t = .true, 
y e ar s _f or _ini t i a 1 _d t = 0.001 

&controls 

okay_to_reduce_gradT_excess = .false. 

use_Type2_opacities = .false. 
use_simple_es_for_kap = .true. 

x-ctrl (l)=ld7 
x_ctrl (2) = 35. 

X -Ctrl (3) = 4. 

X -Ctrl (4) = 0.1 
X-Ctrl (5)=10. 

mixing-length.alpha = 1.89 

min_T_for_acceleration_limited_conv_velocity = 0 

varcontrol-target = le-4 

Dutch-scaling-factor = 0.0 

use_ODE_var_eqn_pairing = . true . 
use_dvdt_form_of_momentum_eqn = . true . 
use_dPrad_dm_form_of_T_gradient_eqn = . true . 

use_compression_outer_BC = . true . 
which.atm.option = ‘ simple.photosphere ’ 
use_zero_dLdm_outer_BC = . true . 

use.artificial-viscosity = .true, 
shock-sp read-linear = 0. 
shock-spread-quadratic = 0.01 

mesh-delta-coeff = 0.8 
min-dq = ld-10 
max-center-cell-dq = 5d-6 
max-surface-cell-dq = ld-10 
log-tau-function-weight = 100 
log-kap-function-weight = 100 

newton-iterations-limit = 7 
iter-for-resid-tol2 = 4 
to 1-r e s i du al-n or m 1 = ld-8 
tol-max-residual 1 = ld-7 

t iny-C orr-Coe f f-li mi t = 999999 

newton_itermin_until_reduce_min_corr_coeff = 999999 
max_age = 350. 
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The choice of max_age is arbitrary. The models were stopped 
roughly after the times shown in Figures [8|10| 

For the model with higher opacity noted in §4.2| we used 
opacityTactor = 100 in addition to the above flags. 

A3 Models From ^ 

These models were run to He core exhaustion using the identical 
inlists as in the first part of Appendix |Al| (the models without heat¬ 
ing) and Dutch_scaling_factor = 2.5, 1.9, 1.75, 1.7, and 0.8 for the 
models in Figure^^fleft to right). Note that decreasing or increas¬ 
ing varcontrolTarget or changing the overshoot/mixing parameters 
can slightly change the time to He core exhaustion and thus the 
stellar mass at that time (the latter because of the effects of stellar 
winds). This can change whether the model is a BSG, RSG, or WR 
star in Figure [TT] However, models with different parameters but 
similar radii at He core exhaustion have similar thermal profiles so 
the thermal profiles in Figure [TT] are more robust than the specific 
model parameters that produce them. 

A4 Extra Energy Deposition 

Our model for extra energy deposition uses the extra_controls hook 
to call a subroutine that deposits a time-independent source of heat¬ 
ing via other_energy as follows: we deposit an energy per unit 
time x_ctrl(l)*Lo as a Gaussian centered at radius x_ctrl(2)*Ro 
with dispersion x_ctrl(3)*Ro- The Gaussian is cutoff below radius 
x_ctrl(4)*R0 (to avoid problems if there is energy deposition too 
deep in the core in some calculations). Finally, we multiply the 
heating by tanh(star_age*x_ctrl(5)) to ensure that a sudden onset 
of super-Eddington heating does not cause numerical problems. 
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